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Introduction 

The  broad,  long-term  objective  of  this  IDEA  proposal  is  to  achieve  optimized  image 
quality  for  DM  within  acceptable  limits  of  radiation  exposure  by  developing  innovative 
approaches  for  controlling  DM  exposures.  These  approaches  entail  using  the  digital 
detector  and  an  artificial  neural  network  to  control  mammographic  exposures.  This 
project's  specific  aims  are  (1)  to  use  short,  low  dose  pre-exposures  of  the  breast  to  create 
“intelligent”  regions  of  interest  that  determine  the  exposure  parameters  for  the  fully 
exposed  image;  and  (2)  to  use  an  artificial  neural  network  to  select  exposure  parameters 
(mAs,  kVp,  and  beam  filtration)  based  on  “intelligent  decisions”  that  optimize  signal-to- 
noise  (SNR)  as  a  function  of  mean  glandular  dose. 


Body:  Progress  on  Research  Tasks  &  Key  Research  Accomplishments 

Task  1.  Determine  optimal  pixel  binning  factor  and  size  of  exposure-controlling 
ROI. 

The  ideal  exposure-controlling  ROI  (ECROI)  is  large  enough  so  that  small,  high-density 
structures  such  as  calcifications  are  not  selected  to  determine  exposure  parameters,  but 
small  enough  so  that  the  image  is  not  smoothed  excessively.  To  determine  an  optimal 
ECROI  size,  we  surveyed  FFDM  images  from  our  clinical  image  database,  evaluating  3- 
D  surface  plots  to  determine  where  the  “peaks”  of  low  exposure  areas  are  located  and  the 
range  of  the  super-pixel  values  in  the  peaks.  This  provided  information  on  the  effects  of 
the  location  and  size  of  the  ECROI.  We  evaluated  the  effect  of  the  ECROI  size  on  the 
calculation  of  the  SNR  (Figures  1-7).  Next,  we  developed  software  to  analyze  selected 
ROIs  (super-pixels)  of  digitally  acquired  mammograms.  The  primary  two  questions  to  be 
addressed  by  these  analysis  tools  are: 

1)  Which  area  of  the  digitally  obtained  mammogram  contains  the  area  of  greatest 
radiographic  attenuation  (the  ECROI)? 

2)  What  portion  of  the  pixel  variance  in  the  ECROI  is  a  result  of  differential 
breast  attenuation  (signal)  as  opposed  to  random  fluctuations  (noise)? 

With  respect  to  question  1),  the  immediate  objectives  of  the  software  programs  are: 

a)  Evaluate  the  correlation  between  the  number  of  regions  of  interest  (ROIs) 
sampled  and  the  correct  identification  of  the  ROI  containing  the  area  of  greatest 
radiographic  attenuation  (the  ECROI) 

b)  Evaluate  the  correlation  between  ROI  size  and  correct  identification  of  the 
ECROI 

c)  Examine  the  effect  of  pixel  binning  on  correct  identification  of  the  ECROI 

The  selection  of  the  pixel  binning  factor  involves  tradeoffs  between  readout  speed  and 
information  loss.  We  evaluated  binning  factors  of  8  (to  form  320p  x  320p  super-pixels) 
and  16  (giving  640p  x  640p  super-pixels)  using  images  of  the  test  phantom  acquired  at 
low  exposures.  Three-dimensional  pixel-value  surface  plots  were  compared  to  similarly 
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binned  plots  made  using  full  exposures  (Figures  1-4).  The  location  of  areas  having  low 
ADU  values  was  the  same  in  the  low-  and  high-exposure  plots,  indicating  that  these 
binning  factors  did  not  result  in  unacceptable  image  smoothing.  We  then  examined  ROIs 
of  800  x  800  pixels  (3200p  X  3200p  pixels)  (Figures  5-7)  and  1600  x  1600  pixels 
(64,000p  X  64,000p  super-pixels).  For  the  3200p  X  3200p  pixel  ROI,  a  6  x  8  grid  of 
ROIs  was  established.  In  our  experiments  (Figures  5  &  6),  the  ROI  positions  were  held 
fixed  while  their  side  lengths  varied  between  minimum  values  of  approximately  1 .7%  of 
the  short  dimension  of  the  image,  to  their  maximum  value,  equal  to  the  center-to-center 
ROI  spacing.  Average  pixel  value  in  the  48  ROIs  (Figure  5)  were  plotted  as  a  function  of 
ROI  side  length,  in  units  of  binned  pixels  (Figure  6).  We  then  performed  similar 
experiments  for  64K  p  X  64K  p  ROIs.  These  data  were  then  utilized  for  task  2 
experiments  (below)  to  calculate  the  contribution  of  signal  to  the  total  variance  of  the 
ECROIs. 

Task  2.  Identify  the  contribution  of  signal  to  the  total  variance  in  the  ECROI. 

This  entails  distinguishing  between  contributions  to  the  variance  within  the  ECROI  from 
breast  structure  and  from  noise  (e.g.,  x-ray  quantum  noise  or  aquantum,  and  additive/system 
noise  or  cre)  in  the  pre-exposure  images.  We  have  used  two  pre-exposures  to  permit  pixel 
fluctuations  arising  from  the  signal  relating  to  breast  structure  to  be  distinguished  from 
those  arising  from  noise.  The  ECROI  is  determined  from  the  first  pre-exposure.  Only 
detector  modules  containing  the  ECROI  are  read  out  on  the  second  pre-exposure.  Pre¬ 
exposure  times  and  initial  kVp  are  scaled  according  to  compressed  breast  thickness  (>6.5 
cm:30  kVp;  6  cm:29  kVp;  5.5  cm:28  kVp;  5  cm:27  kVp;  4.5  cm:26  kVp;  ^4.0  cm:25 
kVp).  Initial  choices  for  pre-exposure  times  and  kVp  are  based  on  values  calculated  to 
result  in  a  MGD  of  approximately  150  mrads  to  a  50%  fat  /  50%  fibroglandular  breast  in 
an  exposure  time  of  approximately  1  sec. 

Criteria  for  optimization  of  tube  voltage  and  external  filtration  in  digital  mammography 
differ  from  those  used  in  screen-film  mammography.  This  is  because  the  separation  of  the 
processes  of  acquisition  and  display  in  the  former  permits  contrast  of  individual 
structures  to  be  adjusted  when  the  image  is  displayed.  It  is  therefore  possible  to  detect 
objects  with  low  subject  contrast  provided  that  the  image  signal  to  noise  ratio  (SNR)  is 
adequate.  Thus,  rather  than  maximization  of  contrast  within  the  constraints  of  acceptable 
film  darkening  and  patient  dose,  beam  optimization  in  digital  mammography  requires 
maximization  of  the  image  SNR,  constrained  by  acceptable  patient  dose. 

To  identify  optimum  technique  factors,  the  following  figure  of  merit  was  used: 

FOM  =  (SNR)2/MGD, 

Where  MGD  is  mean  dose  to  the  glandular  portion  of  the  breast,  and  the  SNR  is  Rose’s 
suggested  minimum  value  of  5  [1]. 

Using  simple  phantoms  with  known  signals  (disks,  holes,  etc.)  that  result  in  signal- 
induced  variances  that  are  easily  determined,  we  calculated  the  ECROI  SNR.  Breast 
tissue  equivalent  material  corresponding  to  3  different  fat/glandular  ratios  was  used  to 
simulate  a  range  of  breast  thicknesses  and  densities.  Look-up  tables  were  generated  with 
starting  pre-exposure  time  and  kVp  for  varying  breast  composition/thickness  (288  test 
conditions)  as  input  data.  These  data  were  further  refined  in  task  3  experiments  using 
5 


Annual  Report:  t_ . . 

L.L. Fajardo,  M.D. 

anthropomorphic  breast  phantoms  that  more  closely  approximate  the  spectral  content  of 
human  breast  tissue  (see  below). 


Task  3.  Determine  how  breast  structure  relates  to  measured  variances.  We  selected 
a  broad  range  of  breast  thickness  and  radiodensity  from  our  Digital  mammography 
patient  image  database  to  empirically  determine  the  contribution  to  the  total  variance  that 
arises  from  breast  structure  (signal).  In  order  to  obtain  the  best  possible  estimate  of  the 
variance  associated  with  normal  breast  parenchyma,  we  selected  ROIs  from  sub-areas  of 
mammograms  with  obviously  high  SNR.  In  those  ROIs,  the  average  pixel  value  and 
system  gain  was  used  to  calculate  the  x-ray  quantum  variance,  aquantum .  Using  the 

measured  system  noise,  ae,  the  contribution  to  the  total  variance,  a  total,  from  breast 
tissue  was  calculated  from:. 


asignal  ~  CTtotal  CTquantum  CTe 

After  the  2  pre-pulses,  we  calculated  cr  noise  and  a  signal  and  used  the  empirically 
determined  relationship  between  a2Signai  and  Imax  -  Imin  to  infer  the  magnitude  of  the 
signal.  The  ratio  of  the  inferred  signal  to  onoise  is  the  SNR  following  the  pre-exposure. 
This  was  then  used  to  select  exposure  parameters  for  the  full  exposure.  If  only  the 
exposure  time  is  to  be  selected  (phototiming  mode),  this  is  straightforward  since  Imax  - 
Imin  is  directly  proportional  to  mAs  and  since  the  quantum  noise,  oq,  increases  as  the 

square  root  of  the  exposure  time.  In  the  CCD-based  TREX  system,  the  system  noise  (ae) 
increases  as  the  square  root  of  the  exposure  time  because  of  the  accumulation  of  thermal 
electrons  in  the  CCD  pixels.  This  noise  was  previously  characterized  over  a  large  range 
of  exposure  times. 


Results  of  Phantom  Experiments  Performed  for  Tasks  2  &  3: 

Following  are  example  plots  of  measured  signal,  SNR,  and  the  FOM,  as  a  function  of 
kVp,  for  a  70/30  fibroglandular/fat  equivalent  phantom  composition,  using  the  signal 
from  the  microcalcification-equivalent  material.  These  data  were  obtained  from  the  DM 
system  at  Johns  Hopkins.  The  results  of  these  studies  were  presented  at  the  5th 
International  Workshop  on  Digital  Mammography,  Toronto,  CA,  June  11-14, 2000. 
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Figure  1 


Figure  2 


Normalized  FOM,  70/30,  Calcs  (Mo  filter) 


20  25  30  35  40 


kVp 

Figure  3 


Nine  different  phantoms  were  assembled  to  simulate  breasts  of  three  different  thicknesses 
(3  cm,  5  cm,  and  7  cm),  and  three  different  attenuation  equivalent  adipose/fibroglandular 
mass  ratios  (30/70,  50/50,  and  70/30).  All  blocks  of  a  given  phantom  had  the  same 
adipose  /fibroglandular  ratio,  except  for  two  5  mm  thick  blocks,  common  to  all  phantoms, 
that  are  100%  adipose  equivalent.  These  blocks  were  placed  at  the  top  and  bottom  of  the 
stack  to  simulate  skin  (see  figure  4).  In  each  phantom  stack  assembled,  the  centrally 
located  block  in  the  stack  (the  signal  block)  contained  a  series  of  test  objects.  For  the  data 
reported  here,  the  test  objects  of  interest  were  two  stepwedges,  one  each  of  calcification 
equivalent  and  mass  equivalent  material.  The  mass  equivalent  stepwedge  has  the  same  x- 
ray  attenuation  as  100%  glandular  equivalent  material,  and  the  microcalcification 
equivalent  step  wedge  is  composed  of  calcium  carbonate.  Figure  5  is  a  schematic  of  a 
signal  block  showing  the  dimensions  of  the  block  and  step  wedges  (other  test  objects 
present  in  the  signal  block  have  been  omitted  for  clarity).  The  thickness  of  all  signal 
blocks  is  2  cm.  Images  were  obtained  in  manual  mode  with  the  phantoms  positioned  at 
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the  chest  wall  edge  of  the  receptor,  centered  left  to  right.  Exposure  time  was  selected  to 
give  approximately  the  same  average  pixel  value  in  the  phantom  background  area  for 
each  phantom/technique  combination.  For  each  combination  two  images  were  obtained 
with  identical  exposure  times  for  the  purpose  of  image  subtraction,  taking  care  not  to 
move  the  phantom  between  the  two  exposures.  At  each  site,  entrance  exposures 
(mR/mAs)  and  half  value  layers  (HVLs)  were  measured  for  each  target/filter/kVp 
combination  used. 


Figure  4:  Side  view  of  a  5  cm  thick  phantom,  comprised  of  one  2  cm  thick  signal  block, 
two  1  cm  thick  blank  blocks,  and  two  0.5  cm  thick  skins. 
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Figure  5:  Schematic  diagram  of  a  signal  block 


ECROI  Analysis: 

Signal  was  defined  as  the  difference  between  the  average  pixel  values  in  a  region  of 
interest  (ROI)  centered  on  an  individual  step  (but  not  including  the  step  boundaries),  and 
an  equal  sized  ROI  located  immediately  adjacent  to  the  step,  but  containing  only 
background.  To  quantify  the  image  noise,  the  two  images  of  a  given  phantom,  obtained  at 
a  common  technique,  were  subtracted.  Image  subtraction  was  performed  to  remove  fixed 
pattern  noise  associated  with  phantom  defects,  detector  non-uniformity,  and  the  heel 
effect.  Noise  in  a  single  image  was  defined  as  the  rms  pixel-to-pixel  fluctuations  in  an 
ROI  of  1 1 09  x  5 1 1  pixels  in  the  difference  image,  divided  by  the  square  root  of  two. 
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Evaluation  of  acceptable  mean  glandular  dose: 

The  MGD  for  each  phantom  was  calculated  using  its  known  thickness  and  composition, 
and  the  measured  HVL  and  mR/mAs  values  from  each  DM  system.  For  Mo/Mo  and 
Mo/Rh  spectra,  the  parameterized  dose  tables  of  Sobol  and  Wu  were  utilized  to  obtain  the 
glandular  dose  per  unit  exposure  [2].  For  the  W/Al  spectra,  normalized  (to  entrance 
exposure)  MGD  values  were  obtained  from  the  data  of  Stanton  et  al.  [3].  Their  data  were 
extrapolated  to  3  cm  breast  thickness,  and  interpolation  between  their  published  HVL 
curves  was  used  to  obtain  correction  factors  for  the  particular  glandular  volume  fractions 
(0.22, 0.40,  and  0.61,  corresponding  to  glandular  mass  fractions  of  0.30,  0.50,  and  0.70, 
respectively)  used  in  our  study.  For  the  W/Rh  spectra,  the  calculations  of  Boone  were 
utilized,  interpolating  between  his  published  HVL  and  adipose/fibroglandular 
composition  values  [4].  All  FOM  values  were  obtained  by  dividing  the  square  of  the  SNR 
by  the  MGD,  expressed  in  units  of  10'5  Gy  (1  mrad). 

The  measured  HVL  values  for  the  seven  specific  target/filter  combinations  tested  at  the 
three  sites,  as  a  function  of  kVp,  are  shown  in  Figure  3.  Figure  4  shows  the  corresponding 
normalized  MGD,  DgN,  calculated  for  each  of  the  seven  spectra,  plotted  versus  the 
measured  HVL.  Similarly,  Figure  8  shows  DgN  for  each  target/filter  combination  tested, 
plotted  versus  kVp.  The  general  tradeoff  between  loss  of  contrast  and  reduction  in  MGD 
with  increasing  kVp  is  illustrated  in  Figure  9.  In  this  example,  the  measured  contrast  of 
the  0.3  mm  thick  (thickest)  calcification  step  in  shown  for  the  5  cm  thick,  50/50  phantom. 

SNR  versus  kVp,  and  the  corresponding  FOM  values  vs.  kVp  have  been  determined. 
Figures  10  and  1 1  show  the  results  obtained  for  the  300  micron  thick  step  of  the 
calcification  stepwedge  in  the  three  50/50  composition  phantoms.  To  illustrate  the 
applicability  of  these  data  to  objects,  the  dependence  of  the  FOM  on  the  step  thickness 
for  both  types  of  stepwedges  is  presented  in  Figures  12  and  13.  These  data  are  from 
images  obtained  using  a  Mo/Mo  target/filter  combination  to  image  a  5  cm  thick,  50/50 
composition  phantom.  Finally,  Figures  14-16  illustrate  the  effect  on  the  FOM  of  changing 
breast  composition,  holding  breast  thickness  fixed.  Signals  were  calculated  using  the  10 
mm  thick  mass  equivalent  step.  , - 

1  1  HVLvskVp 
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Figure  7:  Normalized  mean  glandular 
dose  versus  HVL,  assuming  a  5  cm  thick, 
50/50phantom 


Dose  and  Contrast  vs  kVp 
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Figure  8:  Normalized  mean  glandular 
dose  vs.  kVp  assuming  a  5  cm  thick, 
50/50  phantom 


Figure  9:  Dose  and  contrast  versus  kVp 
using  the  0.3  mm  calcification  step  in  a  5 
cm  thick,  50/50  phantom 
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SNR  vs  kVp:  50/50 


kVp 


Figure  10:  SNR  vs.  kVp. 


FOM  vs  kVp  :  50/50 


Figure  1 1 :  FOM  vs.  kVp. 


Figure  12:  FOM  values  for  the  five  steps  of  the  mass 
stepwedge,  normalized  by  the  average  value  for  each 
step.  The  average  FOM  values  ranged  from  0.2  (step 
0)  to  0.01  l(step  4).  Imaging  data  from  the  5  cm  50/50 
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Figure  13:  FOM  values  for  the  four  steps 
of  the  calcification  stepwedge, 
normalized  by  the  average  value  for  each 
step.  The  average  FOM  values  ranged 
from  1 .4  (step  0)  to  0.64  (step  3).  Data 
are  imaging  the  5  cm  50/50  phantom. 


Figure  14:  FOM  vs  kVp  for  3  cm  thick 
phantoms  of  three  compositions. 
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Figure  15:  FOM  vs  kVp  for  5  cm  thick 
phantoms  of  three  compositions 
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Figure  16:  FOM  vs  kVp  for  7  cm  thick 
phantoms  of  three  compositions. 


Task  4.  Optimize  mAs,  kVp,  and  Filter  Selection  Based  on  the  Results  of  the  Pre¬ 
exposures. 

Optimization  is  based  on  maximization  of  the  figure  of  merit:  FOM=(SNR2/MGD).  We 
have  used  the  added  information  provided  by  the  spatially  varying  signal  levels  in  the 
pre-exposure  ROIs  to  calculate  a  more  refined  estimate  of  the  MGD  than  is  possible 
using  simply  compressed  thickness.  To  do  this,  the  mean  pixel  value  in  each  ROI  is  used 
to  determine  the  relative  transmission  through  the  breast  at  that  location.  In  regions  of 
uniform  thickness  (all  regions  whose  entrance  surface  is  in  direct  contact  with  the 
compression  paddle),  variations  in  x-ray  transmission  are  due  to  variations  in  breast 
composition.  The  region  of  the  breast  above  each  ROI  is  assumed  to  consist  of  two  skin 
and  subcutaneous  layers  with  a  uniform  mixture  of  adipose  and  fibro-glandular  tissue  in 
between.  The  x-ray  transmission  through  each  region  determines  an  adipose/ 
fibroglandular  composition  ratio.  Thus,  for  a  given  kVp,  filter,  and  mAs,  the  MGD  for 
the  region  of  the  breast  above  each  ROI  can  be  calculated  separately,  and  the  result 
summed  to  obtain  an  average  MGD. 

Preliminary  Results:  Task  4: 

Software  has  been  developed  for  analysis  of  digitally  acquired  mammograms.  The 
primary  two  questions  to  be  addressed  by  these  analysis  tools  may  be  generalized  as 
follows: 

1)  Which  area  of  the  digitally  obtained  mammogram  contains  the  area  of  greatest 
radiographic  attenuation  (the  ECROI)? 

2)  What  portion  of  the  pixel  variance  in  the  ECROI  is  a  result  of  differential  breast 
attenuation  (signal)  as  opposed  to  random  fluctuations  (noise)? 

With  respect  to  question  1),  the  immediate  objectives  of  these  programs  are: 

1)  Evaluate  the  correlation  between  the  number  of  regions  of  interest  (ROIs) 
sampled  and  the  correct  identification  of  the  ROI  containing  the  area  of  greatest 
radiographic  attenuation  (the  ECROI) 

2)  Evaluate  the  correlation  between  ROI  size  and  correct  identification  of  the 
ECROI 


13 


Annual  Report: 


L.LFajardo,  M.D. 


3)  Examine  the  effect  of  pixel  binning  on  correct  identification  of  the  ECROI 

Figure  17  below  shows  an  example  of  a  digitally  obtained  mammogram  in  which  a  6  x  8 
grid  of  ROIs  has  been  established.  The  figure  shows  the  mammogram,  labeled  with 
numbers  whose  locations  correspond  to  the  centers  of  the  ROIs. 


Figure  17:  A  digitally  acquired  mammogram, 


annotated  with  the  locations  of  a  6  x  8  square  grid 
of  ROIs.  The  numerals  are  located  at  the  centers 
of  the  square  ROIs.  The  ROI  positions  are  held 
fixed,  while  their  side  lengths  are  varied  between 
minimum  values  of  approximately  1 .7%  of  the 
short  dimension  of  the  image,  to  their  maximum 
value,  equal  to  the  center-to-center  ROI  spacing. 


Figure  18:  Average  pixel  value  in  the  48  ROIs  of  Figure 
1,  plotted  as  a  function  of  ROI  side  length,  in  units  of 
binned  pixels.  For  the  example  shown,  the  ECROI  (that 
with  the  lowest  average  pixel  value)  is  the  same  for  all 
ROI  sizes  tested.  Its  center  is  at  the  position  labeled  ‘19’ 
in  Figure  17.  ROIs  whose  average  pixel  value  is  always 
large  and  constant  (top  of  the  graph)  are  located  outside 
the  region  of  the  breast. 


Figure  18  is  a  plot  of  the  mean  values  of  each  ROI,  plotted  as  a  function  of  ROI  size. 
The  ROI  with  the  lowest  value  is  the  ECROI.  For  the  example  shown,  the  ECROI  is 
ROI  #19.  The  contour  plot  of  the  pixel  values,  shown  below  in  Figure  19,  verifies 
that  this  is  in  fact  the  most  highly  attenuating  portion  of  the  breast. 
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Figure  19:  Contour  plot  of  the  pixel  values  in  the  digital 
mammogram  of  Figure  1. 


Ongoing  work 

We  are  now  in  the  process  of  running  our  analysis  tools  using  digital  mammograms 
obtained  from  breasts  with  a  variety  of  compositions  and  sizes.  We  will  evaluate  the  effect 
of  radiographic  heterogeneity  on  the  selection  of  the  ECROI,  as  well  as  possible  spurious 
results  from  microcalcification  clusters  or  cysts. 

We  are  also  developing  analysis  programs  for  determination  of  the  pixel  variance  due  to 
differential  breast  attenuation,  and  are  characterizing  intrinsic  detector  noise. 

Using  the  data  obtained  from  our  experiments  to  date,  we  have  implemented  an  automatic 
exposure  control  device  onto  the  digital  mammography  system  that  will  be  refined  during 
the  remaining  6  months  of  support  from  this  grant.  Table  1  illustrates  the  technique  factors 
and  dose  generated  by  the  AEC  of  a  conventional  screen-film  unit  compared  with  those 
generated  by  the  AEC  implemented  on  our  digital  mammography  system  for  patients 
undergoing  both  studies.  Similar  or  reduced  mean  glandular  doses  were  produced  by  the 
AEC  implemented  on  the  digital  system,  with  one  exception. 
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Reportable  Outcomes: 

Abstracts  &  Presentations: 

The  1999  Radiological  Society  of  North  America,  Chicago,  IL,  November  27-  December  3, 
1999:  “Development  of  a  quality  control  system  for  full-field  digital  mammography”. 

MJ  Yaffe,  MB  Williams,  LT  Niklason,  GE  Mawdsley,  AD  Maidment,  Radiology  209(P)  (1999) 
160. 

The  5th  International  Workshop  on  Digital  Mammography,  Toronto,  CA,  June  11-14,  2000: 
“Beam  Optimization  for  Digital  Mammography”.  MB  Williams,  M  More,  V  Venkatakrishnan,  L 
Niklason,  MJ  Yaffe,  G  Mawdsley,  A  Bloomquist,  A  Maidment,  D  Chakraborty,  C  Kimme- 
Smith,  LL  Fajardo. 

Publication: 

Williams  MB,  Mangiafico  PA,  Simoni  PU.  Noise  power  spectra  of  images  from  digital 
mammography  detectors.  Med  Phys  2000:6(7);1279-1293. 


Conclusions: 

The  analysis  of  SNR  and  FOM  as  a  function  of  kVp,  shown  in  Figures  10  and  11,  indicates  that 
although  the  image  SNR  tends  to  decrease  monotonically  for  all  systems  with  increasing  kVp, 
the  accompanying  MGD  reduction  results  in  fairly  flat  FOM  curves.  This  is  primarily  due  to  tube 
loading,  since  it  was  not  possible  to  obtain  the  same  exit  exposure  at  all  kVps  (that  is,  the  tube 
output  was  insufficient  to  compensate  for  the  lower  transmission  through  the  phantoms).  Thus 
the  falling  SNR  (and  the  falling  MGD)  with  decreasing  kVp  are  really  consequences  of  falling 
exposure. 

For  a  given  phantom/technique  combination,  the  SNR,  and  thus  the  magnitude  of  the  FOM, 
increases  with  increasing  step  thickness  for  both  types  of  stepwedge.  However,  the  shape  of  the 
FOM  vs.  kVp  curves  for  a  given  target/filter/phantom  combination  are  essentially  independent  of 
step  thickness,  and  are  similar  for  mass  and  calcification  equivalent  signals.  This  is  illustrated  by 
the  example  shown  in  Figures  12  and  13.  This  implies  that  the  result  of  the  optimization  is  not 
sensitively  dependent  on  signal  amplitude. 

Figures  14-16  illustrate  a  clear  advantage  to  using  rhodium  filtration  for  thin  breasts,  but  that  for 
breasts  5  cm  or  thicker,  aluminum  filtration  becomes  increasingly  advantageous.  Similar 
statements  can  be  made  for  the  molybdenum  target  systems  tested,  where  molybdenum  filtration 
was  superior  for  3  cm  phantoms  of  all  compositions,  but  rhodium  filtration  produced  better 
results  for  5  and  7  cm  thick  phantoms  of  all  compositions.  These  data  suggest  that  the  choice  of 
external  filtration  is  potentially  more  significant  in  determination  of  the  overall  FOM  of  a  DM 
system  than  is  choice  of  tube  voltage. 

Fahrig  and  Yaffe  developed  a  model  for  optimizing  spectral  shape  in  digital  mammography,  and 
used  it  to  calculate  kVp  values  producing  maximum  SNR  at  a  fixed  dose  for  W  and  Mo  spectra 
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[5].  They  found  that,  for  a  fixed  MGD  of  0.6  mGy  (60  mrads),  the  peak  SNR  occurred  in  the  24- 
31  kVp  range  (W  spectrum)  and  25-29  kVp  range  (Mo  spectrum)  for  4  -  8  cm  breast  thickness, 
and  50/50  breast  composition.  Their  results  were  the  same,  whether  the  lesion  type  modeled  was 
infiltrating  ductal  carcinoma  or  microcalcification. 

Jennings  et  al.  used  a  computational  approach  to  identify  maximum  FOM  values  (FOM  = 
SNR2/MGD)  for  a  variety  of  target/filter  combinations,  and  breast  thicknesses  [6].  They  found 
that  for  a  Mo/Mo  beam  used  to  image  3-6  cm,  50/50  breasts,  the  FOM  peaks  at  27-28  kVp,  and 
changes  slowly  with  changing  kVp  near  the  peak  values.  Very  similar  FOM  vs.  kVp  curves  were 
obtained  for  Mo/Mo,  Mo/Rh,  and  W/Al  spectra,  applied  to  6  cm  thick,  50/50  composition 
breasts.  The  general  trends  in  our  data  appear  to  be  consistent  with  those  of  these  previous 
studies. 

Implementation  of  the  first  generation  of  automatic  exposure  control  device  on  the  digital 
mammography  system  results  in  images  having  satisfactory  image  quality,  SNR  and  mean 
glandular  dose.  For  16  patients  undergoing  both  conventional  mammography  and  digital 
mammography,  the  implemented  AEC  for  the  digital  system  performed  satisfactorily. 
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ABSTRACT 

Because  image  noise  constitutes  the  fundamental  limitation  of  perceptibility  in  a  digital  radiograph,  noise  characterization 
through  estimation  of  the  noise  power  spectrum  (NPS)  is  a  central  component  of  the  evaluation  of  digital  x-ray  systems. 
We  begin  with  a  brief  review  of  the  fundamentals  of  NPS  theory  and  measurement,  derive  explicit  expressions  for 
calculation  of  the  one-  and  two-dimensional  (1-d  and  2-d)  NPS,  and  discuss  some  of  the  considerations  and  tradeoffs 
when  these  concepts  are  applied  to  digital  systems.  Measurements  of  the  NPS  of  two  detectors  for  digital  mammography 
are  presented  to  illustrate  some  of  the  implications  of  the  choices  available.  One  mammographic  detector  is  part  of  an 
early  commercial  system,  currently  under  initial  clinical  evaluation.  The  other  detector  is  being  developed  under  a 
collaborative  effort  between  Brandeis  University  (Waltham,  MA)  and  the  University  of  Virginia  (Charlottesville,  VA).  Two 
dimensional  (2-d)  noise  power  spectra  obtained  over  a  range  of  input  fluences  exhibit  pronounced  asymmetry  between 
the  orthogonal  frequency  dimensions.  The  2-d  spectra  of  both  systems  also  demonstrate  dominant  periodic  structures 
both  on  and  off  the  primary  frequency  axes.  Our  measurements  indicate  that  for  these  systems,  1-d  spectra  useful  for  the 
purpose  of  DQE  calculation  may  be  obtained  from  thin  cuts  through  the  central  portion  of  the  calculated  2-d  NPS.  On  the 
other  hand,  low  frequency  spectral  values  do  not  converge  to  an  asymptotic  value  with  increasing  slit  length  when  1-d 
spectra  are  generated  using  the  scanned  synthetic  slit  method. 


I.  INTRODUCTION 

Because  image  contrast  can  be  manipulated  during  the  display  of  digitally  acquired  radiographic  images,  noise  imposes 
the  fundamental  limitation  to  object  perceptibility  in  a  digital  radiograph.  Thus  noise  characterization  plays  an  increasingly 
central  role  in  the  evaluation  of  system  performance  in  medical  imaging.  Well-developed  methods  used  for  noise 
characterization  in  film-based  systems  are  currently  being  reassessed  and  reformulated  to  accommodate  both  the 
different  nature  of  the  noise  and  to  take  advantage  of  the  analytical  tools  available  with  digital  processing.  An  indication  of 
both  the  current  interest  level  and  lack  of  consensus  regarding  noise  characterization  is  the  recent  establishment  of  an 
AAPM  Task  Group  (#16)  for  Standards  for  Noise  Power  Spectrum  Analysis. 

The  noise  power  spectrum  (NPS)  is  a  spectral  decomposition  of  the  variance.  As  such,  the  NPS  of  a  digital  radiographic 
image  provides  an  estimate  of  the  spatial  frequency  dependence  of  the  pixel-to-pixel  fluctuations  present  in  the  image. 
Such  fluctuations  are  due  to  the  shot  (quantum)  noise  in  the  x-ray  quanta  incident  on  the  detector,  and  any  noise 
introduced  by  the  series  of  conversions  and  transmissions  of  quanta  in  the  cascaded  stages  between  detector  input  and 
output.  Examples  of  the  latter  are  gain  variances  in  the  conversion  of  x-ray  quanta  to  light  quanta  in  a  phosphor  (or  to 
electron-hole  pairs  in  a  solid  state  detector),  statistical  fluctuations  in  the  transmission  of  optical  quanta  between  a 
scintillator  and  a  photodetector,  and  additive  noise  sources  such  as  preamplifier  noise.  Because  it  gives  information  on 
the  distribution  in  frequency  space  of  the  noise  power,  the  NPS  is  a  much  more  complete  description  of  image  noise  than 
is  quantification  of  integrated  (total)  noise  via  simple  measurement  of  the  rms  pixel  fluctuations.  The  total  variance  can  be 
obtained  by  integrating  the  NPS  over  spatial  frequency. 

In  this  paper  we  discuss  aspects  of  NPS  estimation  as  applied  to  digital  radiographic  systems.  We  begin  with  a  general 
review  of  pertinent  definitions  and  relations,  and  discuss  the  most  important  considerations  when  applying  spectral 
estimation  techniques  to  digital  radiographic  systems.  Mathematical  derivations  sufficient  to  provide  the  reader  with  tools 
for  quantitative  application  of  the  concepts  are  given,  with  more  detailed  development  available  in  the  cited  references. 
These  techniques  are  then  applied  to  the  analysis  of  the  noise  properties  of  two  developmental  CCD-based  detectors  for 
full  field  digital  mammography. 


II.  BACKGROUND 

Linear  systems  analysis  can  be  usefully  applied  to  any  radiographic  system  whose  transfer  characteristic  between  input 
and  output  is  linear  ( or  linearizable)  and  stationary  (the  transfer  properties  are  independent  of  location  at  the  input 
surface)  [H.H.  Barrett  and  W.  Swindell,  Radiological  Imaging:  The  theory  of  image  formation,  detection,  and  processing, 
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’  New  Yorkv  Academic  Press,  1981.,  K.  Rossmann,  “The  spatial  frequency  spectrum:  A  means  for  studying  the  quality  of 
radiographic  imaging  systems",  Radiology  90(1):1-13, 1968.].  In  this  formalism,  signal  and  noise  are  typically  decomposed 
into  their  Fourier  components  (although  other  bases  can  be  used),  and  input  and  output  relations  are  expressed  in  terms 
of  the  frequency  space  Fourier  amplitudes.  Digital  radiographic  systems  such  as  the  one  considered  here  usually  exhibit 
linear  responses  to  the  incident  x-ray  fluence  over  several  orders  of  magnitude.  The  following  analysis  assumes  that  the 
detector  is  operated  within  its  linear  range.  Because  in  a  digital  detector  the  analog  (presampling)  signal  and  noise 
propagated  through  the  detector  are  sampled  by  the  digital  matrix  at  discrete  points,  the  system  response  is  not  strictly 
shift  invariant  (unless  shifts  are  by  an  integer  number  of  pixels).  The  effects  of  digital  sampling  on  the  NPS  are  described 
below. 

A  random  process  such  as  the  noise  fluctuations  in  an  image  is  considered  (wide-sense)  statistically  stationary,  if  its 
autocorrelation  function  (see  below)  is  independent  of  the  particular  data  sample  (e.g.  the  particular  location  in  the  image) 
used  to  obtain  it.  This  is  only  strictly  true  in  a  radiographic  image  if  a)  noise  samples  are  obtained  from  spatially  uniform 
exposures,  and  b)  the  detector  adds  no  spatially  fixed  correlated  noise.  Dobbins  has  shown  that  the  presence  of  a 
deterministic  signal  in  the  image  data  used  to  estimate  the  NPS  has  no  effect  on  the  resultant  spectral  estimates  [J.T. 
Dobbins,  “Effects  of  undersampling  on  the  proper  interpretation  of  modulation  transfer  function,  noise  power  spectra,  and 
noise  equivalent  quanta  of  digital  imaging  systems",  Med.  Phys.  22(2),  171-181  (1995)].  One  other  important  assumption 
concerning  the  noise  is  whether  or  not  it  is  ergodic.  The  noise  in  an  image  is  ergodic  if  expectation  values  (averages) 
obtained  from  data  samples  at  various  locations  in  the  image  are  equivalent  to  ensemble  averages  obtained  from 
repeated  measurements  under  identical  conditions  at  a  single  location.  Ergodicity  implies  stationarity,  but  not  necessarily 
vice  versa.  The  assumption  of  ergodicity  of  the  noise  in  uniform  irradiation  radiographic  images  is  typically  made,  and 
permits  averaging  of  spectra  obtained  from  several  areas  of  the  image. 

In  a  cascaded  system  such  as  a  digital  mammography  detector,  the  output  of  one  stage  constitutes  the  input  to  the 
following  stage.  In  a  linear,  shift  invariant  cascaded  system,  each  of  the  independent  stages  are  also  linear  and  shift 
invariant.  At  the  output  of  any  given  stage  of  a  cascaded  detector  system,  the  minimum  possible  value  of  the  NPS  is  that 
set  by  the  (uncorrelated)  Poisson  variance  of  the  output  quanta,  which  is  equal  to  the  number  of  quanta.  Building  on  the 
work  of  Rossmann  [K.  Rossmann,  “Spatial  fluctuations  of  x-ray  quanta  and  the  recording  of  radiographic  motttle”,  Am.  J. 
Roentgenol.  90,  863-869  (1963).]  and  Kemperman  and  Trabka  [J.H.B.  Kemperman  and  E.A.  Trabka,  “Exposure 
fluctuations  in  model  photographic  emulsions.  II.  Statistics  due  to  intensifying  screens",  Photogr.  Sci.  Eng.  25, 108-111 
(1981).],  Rabbani  et  al.  have  developed  expressions  for  the  propagation  of  the  NPS  through  cascaded  linear  systems  [M. 
Rabbani,  R.  Shaw,  and  R.  Van  Metter,  “Detective  quantum  efficiency  of  imaging  systems  with  amplifying  and  scattering 
mechanisms’,  Opt.  Soc.  Am.  A.  4(5),  895-901  (1987).].  Using  moment  generating  functions,  they  show  that  noise 
generated  during  a  stochastic  amplification  process  in  one  stage  (such  as  the  conversion  of  a  single  absorbed  x-ray 
photon  to  many  visible  photons)  is  effectively  passed  as  signal  to  the  subsequent  stage.  Furthermore,  they  show  that  in  a 
stochastic  scattering  process,  while  correlated  components  of  the  NPS  are  reduced  by  the  square  of  the  modulation 
transfer  function  (MTF)  associated  with  the  blurring  process,  uncorrelated  Poisson  noise  is  transferred  independently  of 
the  MTF.  This  stochastic  type  of  blur  has  a  different  effect  on  the  NPS  than  deterministic  blur,  in  which  both  correlated  and 
uncorrelated  components  of  the  NPS  of  the  preceding  stage  are  attenuated  by  the  square  of  the  blur  MTF  [I.A. 
Cunningham,  M.S.  Westmore,  and  A.  Fenster,  “Effect  of  finite  detector  element  width  on  the  spatial  frequency  dependent 
detective  quantum  efficiency",  Proc  SPIE  2432  143-151  (1995).].  Because  the  signal  is  also  modulated  by  the  blur  MTF, 
deterministic  blur  does  not  degrade  the  signal  to  noise  ratio.  Visible  photon  scatter  in  a  phosphor  is  an  example  of 
stochastic  blur;  the  effect  of  the  integrating  aperture  of  the  individual  pixels  in  a  digital  system  is  an  example  of 
deterministic  blur. 


A.  Direct  and  Indirect  Methods  of  NPS  Calculation 

1.  The  NPS  is  defined  as  the  Fourier  transform  of  the  autocorrelation  function,  defined  in  one  dimension  as  [R.B. 
Blackman  and  J.W.  T ukey,  The  Measurement  of  Power  Spectra,  (Dover  Publications,  Inc.,  New  York,  1 958),  Chap. 

1]: 

X/2 

C(x)=  lim  j  l(x+T)l  (T)dt  (1) 

x-»°°-X/2 

where  l(x)  is  the  value  (pixel  vajue  for  a  digital  image)  of  the  one  dimensional  image  at  position  x,  and  l*(x)  is  its  complex 
conjugate.  Since  l(x)  is  real,  1=1*.  The  Fourier  transform  of  C(x)  is  then: 


S(u)=  |C(x)e_2,dxiJdx 
00 
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=  lim  —  J  l(T)e+2niTU  Jl(x+x)e_2nixudx[e  2mujdT 

X— >c»  X  -X/2  _oo  ^  ** 

1  X/2 

=  Jim  -^l(u)  J  l(x)e+2niTUdt  (2) 

X->°°  *  -X/2 

Inspection  of  equations  (1)  and  (2)  show  that: 

S(u)  =  FT{C(x)}  =  Jim  l|l(u)|2  (3) 

X — >oo  A 


where  l(u)  =  FT{l{x)}. 

Thus  the  NPS  may  be  calculated  either  from  the  Fourier  transform  of  the  autocorrelation  function  (the  indirect  method),  or 
from  the  square  of  the  modulus  of  the  Fourier  transform  of  the  data  itself  (the  direct  method).  Note  that  since  the  data  l(x) 
are  real,  that  l(-u)=f  (u).  That  is,  the  NPS  at  a  given  negative  frequency  is  equal  in  magnitude  to  the  NPS  at  the 
corresponding  positive  frequency.  Extending  the  expression  of  equation  (3)  to  two  dimensions: 

S(u,v)  =  FT{C(x,y)}  =  lim  -U(u,v)|2  (4) 

X,Y-»oe  XY 

With  the  advent  of  the  Fast  Fourier  Transform  (FFT)  and  fast  computers,  indirect  calculation  of  the  NPS  via  the 
autocorrelation  function  has  largely  been  replaced  by  the  direct  method.  This  paper  will  discuss  only  the  latter  method. 


B.  One  Dimensional  NPS 

The  detective  quantum  efficiency  (DQE)  is  perhaps  the  best  single  descriptor  of  radiographic  detector  performance.  It 
quantifies  how  well  the  detector  is  able  to  transfer  the  signal  to  noise  ratio  (SNR)  inherent  in  the  x-ray  fluence  at  its  input 
to  the  resulting  image.  For  most  two  dimensional  digital  imaging  detectors,  the  DQE  is  a  function  of  at  least  three 
independent  variables:  its  two  spatial  coordinates,  and  the  x-ray  exposure.  It  is  a  function  of  a  single  spatial  coordinate  for 
the  special  case  of  detectors  with  an  isoplanatic  response.  The  DQE  can  be  written  [R.  Shaw,  The  equivalent  quantum 
efficiency  of  the  photographic  process",  J.  Photogr.  Sc.  1 1 , 1 99-204  (1 963),  J.C.  Dainty  and  R.  Shaw,  Image  Science 
(Academic  Press,  New  York,  1974)  pg.  312.]: 


DQE(u,v,<p) 


(<pg)2  C(u,v)MTF2(u,v) 
S(u,v,<p)SNRfn 


(5) 


where  9  represents  the  magnitude  of  the  input  x-ray  fluence  at  the  detector  surface,  g  relates  changes  in  the  zero 
frequency  output  signal  to  that  at  the  input  (i.e.  it  is  the  system  gain),  C(u,v)  is  the  Fourier  decomposition  of  the  input 
signal,  MTF(u,v)  is  the  detector  modulation  transfer  function,  and  S(u,v,9)  is  the  output  noise  power  spectrum  [J.M. 
Sandrik  and  R.F.  Wagner,  “Absolute  measures  of  image  quality  in  magnification  radiography”,  Med.  Phys.  9(4),  540-549 
(1982).].  For  a  monochromatic  input  x-ray  spectrum,  or  in  the  case  of  a  photon  counting  x-ray  detector,  SNRin  is  given  by 
the  magnitude  of  the  input  fluence,  9. 


In  part  because  illustration  of  DQE(u,v,  9)  requires  a  4-dimensional  medium,  DQEs  of  2-d  imaging  systems  are  typically 
presented  one  spatial  frequency  coordinate  at  a  time  for  clarity.  By  definition: 

MTF(u,v)=l2dFT{  psf(x.y)}  I  (6) 

where  psf(x.y)  is  the  spatial  response  of  the  detector  to  a  point  x-ray  input,  and  2dFT  represents  a  two  dimensional 
Fourier  transform.  The  1-d  MTF  is  typically  measured  using  a  narrow  slit  or  edge  [I.A.  Cunningham  and  B.K.  Reid,  “Signal 
and  noise  in  modulation  transfer  function  determinations  using  the  slit,  wire,  and  edge  techniques”,  Med.  Phys.  19(4), 
1037-1044  (1992).],  and  is: 


MTF(u)  =  FT(lsf(x)}  =  FT{  -£-[  esf(x)  ] } 


(7) 


Since  lsf(x)=  |psf(x,y)dy ,  then 
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|  lsf(x)e  2roxudx  =  J{psf(x,y)dy}e  2nixudxj 


=  |2dFT{psf(x,y)}|v_0  =  MTF(u,0) 


(8) 


Equations  (7)  and  (8)  indicate  that  evaluation  of  FT{lsf(x)}  yields  MTF(u,0),  and  thus  the  appropriate  expression  of  the 
NPS  for  calculation  of  the  single  coordinate  DQE  is  S(u,0).  Similar  arguments  pertain  to  MTF(y)  [B.H.  Hasegawa,  The 
Physics  of  Medical  X-Ray  Imaging,  2nd  ed.  (Medical  Physics  Publishing,  Madison,  Wl,  1991)]. 


Up  to  now,  most  characterizations  of  S(u,0)  (or  S(0,v))  of  either  analog  (film-based)  or  digital  radiographic  systems  have 
been  made  using  the  synthesized  slit  method  [W.  Hillen,  U.  Schiebel,  and  T.  Zaengel,  “Imaging  performance  of  a  digital 
storage  phosphor  system”,  Med.  Phys.  14(5),  744-751  (1987).  H.  Roehrig,  L.  Fajardo,  and  T.  Yu,  “Digital  x-ray  cameras 
for  real  time  steroetactic  breast  needle  biopsy”,  Proc.  SPIE  1896,  213-224  (1993).  A.D.A.  Maidment  and  M.J.  Yaffe, 
“Analysis  of  the  spatial-frequency-dependent  DQE  of  optically  coupled  digital  mammography  detectors”,  Med.  Phys. 

21(6),  721-729  (1994).  W.  Zhao,  I.  Blevis,  S  Germann,  and  J.A.  Rowlands,  “Digital  radiography  using  active  matrix 
readout  of  amorphous  selenium:  Construction  and  evaluation  of  a  prototype  real-time  detector”,  Med.  Phys.  24(12),  1834- 
1843  (1997).].  In  the  original  development  of  the  method  in  the  analog  context,  the  optical  aperture  of  the 
microdensitometer  used  to  measure  the  film  optical  density  formed  the  slit.  The  methodology  has  since  been  applied  to 
digital  systems  by  forming  an  effective  slit  from  rectangular  groups  of  contiguous  pixels.  In  either  case,  a  long,  narrow  slit 
is  scanned  in  steps  in  one  direction  over  the  image.  At  each  step  the  average  intensity  in  the  slit  (average  transmitted  light 
for  a  scanning  microdensitometer,  or  average  pixel  value  for  a  digital  system),  is  recorded.  In  this  way,  a  one  dimensional 
data  series  is  generated.  The  squared  modulus  of  the  Fourier  transform  of  the  series,  after  scaling,  gives  the  measured  1- 
d  NPS.  In  the  limit  of  an  infinitely  long  1-d  data  series  in  the  x-direction,  the  resulting  NPS  is  [R.C.  Jones,  “New  method  of 
describing  and  measuring  the  granularity  of  photographic  materials”,  J.  Opt.  Soc.  Am.  45, 799  (1955).  J.C.  Dainty  and  R. 
Shaw,  Image  Science  (Academic,  London,  1974),  pg.  279.]: 


Sss(u)=  JS(u,v)|T(u,vf  dv  (9) 

— oo 

where  T(u,v)  is  the  optical  transfer  function  (OTF)  of  the  slit,  and  where  S^u)  is  the  1-d  NPS  measured  via  the  scanned 
slit  method.  If  the  slit  is  rectangular,  with  x-dimension  w  and  y-dimension  L,  then  its  point  spread  function  is: 


and  therefore  T(u,v)  is: 


where  the  sine  function  is  defined  as 


T(u,v)  =  FT{h(x,y)}  =  sinc(wu)  sinc(Lv) 


sine  (x)  = 


smrtx 

rtX 


(10) 


(11) 

(12) 


Equation  (9)  becomes: 


Sss(u)  =  sinc2(wu)  jS(u,v)sinc2(Lv)dv 


(13) 


If  the  slit  width  is  chosen  sufficiently  narrow,  then  sine  (wu)  ~  1  within  the  range  of  u  over  which  Sss(u)  is  to  be  evaluated. 

In  any  case,  correction  of  the  measured  NPS  by  this  factor  is  straightforward.  In  a  digital  system,  a  synthesized  slit  of  1  x 
Nl  pixels  constitutes  sampling  the  presampling  NPS  (which  already  includes  the  smoothing  effect  of  the  nonzero  pixel 
aperture),  using  a  line  of  NL  delta  functions,  at  an  interval  in  the  scan  direction  equal  to  the  pixel  width,  Ax.  The  NL 
sampled  values  are  then  averaged.  Thus,  sine  (wu)  ->  sine  (0)  =  1 ,  and  L  =  NL  Ay,  where  Ay  is  the  pixel  size  perpendicular 
to  the  scan  direction. 


As  L-*»,  sinc2(Lv)  becomes  sufficiently  narrow  that  S(u,v)  is  approximately  constant  over  the  range  of  v  (near  v=0)  in 
which  sincz(Lv)  is  appreciably  nonzero  [J.C.  Dainty  and  R.  Shaw,  Image  Science  (Academic,  London,  1974),  pg.  279], 
and 


Sss(u)  =  J S(u,v)sinc2  (Lv)dv  =  S(u,0)  Jsinc2  (Lv)dv  = 


S(u,0) 


Thus  the  desired  NPS  isS(u,0)  =  LSss(u). 


(14) 
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.  Substantial  effort  has  been  devoted  to  understanding  the  criteria  for  "sufficiently  large  L”.  Sandrick  and  Wagner  studied 
the  effect  of  varying  slit  length  on  the  1-d  NPS  of  radiographic  screen-film  systems  [J.M.  Sandrick  and  R.F.  Wagner, 
“Radiographic  screen-film  noise  power  spectrum:  variation  with  microdensitometer  slit  length”,  Appl.  Opt.  20(16),  2795- 
2798  (1981)).  The  found  that  the  magnitude  of  the  low  frequency  NPS  is  underestimated  if  the  slit  length  is  insufficient,  but 
that  the  low  frequency  components  approach  a  limiting  value  as  the  slit  length  is  increased.  The  slit  length  required  for  low 
frequency  NPS  values  within  5%  of  the  plateau  value  was  dependent  on  the  particular  screen-film  combination,  but 
ranged  between  a  value  equal  to  the  length  of  the  scan  and  one  nearly  twice  as  long.  Koedooder  et  al.  developed  an 
expression  to  numerically  approximate  the  integral  in  equation  (13)  [K.  Koedooder,  J.  Strackee,  and  H.W.  Venema,  “A 
new  method  for  microdensitometer  slit  length  correction  of  radiographic  noise  power  spectra”,  Med.  Phys.  13(4),  469-473 
(1 986)].  The  technique  is  applicable  only  to  systems  with  rotational  symmetry. 

An  alternate  to  the  scanned  slit  approach  is  to  extract  a  slice  along  one  of  the  primary  spatial  frequency  axes  from  the  2-d 
NPS.  However,  because  these  slices  contain  noise  power  from  the  orthogonal  dimension  near  and  at  zero  frequency,  low 
frequency  trending  effects  in  the  data  can  dominate  their  magnitude  unless  these  trends  are  first  removed.  For  this 
reason,  the  limited  number  of  studies  employing  this  approach  have  used  “thick”  slices  made  by  extracting  NPS  data  near 
to,  but  not  directly  on,  the  primary  axes  [J.T.  Dobbins  III,  D.L.  Ergun,  L.  Rutz,  D.A.  Hinshaw,  H.  Blume,  and  D.C.  Clark, 
“DQE(f)  of  four  generations  of  computed  radiography  acquisition  devices",  Med.  Phys.  22(10),  1581-1593  (1995),  R.J. 
Jennings,  H.  Jafroude,  R.M.  Gagne,  T.R.  Fewell,  P.W.  Quinn,  D.E.  Steller  Artz,  J.J.  Vucich,  M.T.  Freedman,  and  S.K. 

Mun,  “Storage-phosphor-based  digital  mammography  using  a  low-dose  x-ray  system  optimized  for  screen-film 
mammography",  SPIE  2708, 220-232  (1996).].  As  will  be  demonstrated  below,  in  some  digital  detectors  this  approach  can 
overlook  structures  that  occur  only  at  zero  frequency  in  the  orthogonal  dimension. 


C.  Digital  Noise  Power  Spectrum 

In  the  above  expressions,  the  1-d  and  2-d  NPS  are  written  as  functions  of  a  continuous  variable,  I,  that  describes  the 
point-to-point  fluctuations  in  the  image.  In  a  digital  system,  this  continuous  function  of  the  spatial  coordinates  is  sampled 
at  regular  intervals,  Ax.  Thus  I  is  evaluated  at  a  set  of  discrete  locations,  x=nAx,  n=0,1,2,...N,  where  X=NAx. 
Correspondingly,  functions  of  the  spatial  frequency  variable  u  are  now  discrete;  u=kAu,  k=0,1 ,2,...K,  where  KAu  is  the 
maximum  spatial  frequency  sampled.  This  maximum  measurable  frequency,  the  Nyquist  frequency,  is  Un  =  KAu  =  1/(2Ax). 
To  calculate  the  digital  NPS,  the  above  continuous  expressions  must  be  rewritten  in  terms  of  discrete  variables.  With  the 
preceding  definitions,  equation  (3)  can  be  written: 


S(kAu)  =  lim 


AxiVax^-2™^^ 

n=0 


NAx-*ooNAx| 

'KnAx)e-2,ri(kAu)(nAx) 

) 

where  N  is  the  number  of  samples  in  the  interval  X.  Similarly,  equation  (4)  for  the  2-d  NPS  becomes: 


,.  Ax 
=  lim  — 

NAx— N 


(15) 


S(kAu,IAv)  =  lim 
NAx->°° 
MAy-»oo 


Ax  Ay 

N  M 


I  I lKnAx,mAy)e-2*if(kAu)(nAx)+(IAv)(mAy)l 
n=0  m-0 


(16) 


The  experimentally  determined  estimate  of  the  digital  NPS  differs  from  the  underlying  analog  spectrum  in  the  following 
ways: 


a)  the  noise  has  been  convoluted  with  the  rect  (boxcar)  function  given  by  the  pixel  aperture  (deterministic  blurring) 

b)  the  noise  data  has  been  sampled  at  discrete  intervals  given  by  the  pixel  center-to-center  spacing  (equivalent  to 
multiplication  by  a  comb  function) 

c)  the  record  (image  data)  used  to  compute  the  estimate  has  been  truncated  from  its  original  infinite  length  to  a  length 
given  in  each  dimension  by  the  number  of  points  times  the  sampling  interval  (equivalent  to  multiplication  by  a  2-d 
boxcar  function  of  dimensions  NAx  x  MAy) 


The  result  of  a)  is  the  multiplication  of  the  analog  NPS  with  the  absolute  square  of  the  2-d  sine  function  that  is  the  Fourier 
transform  of  the  pixel  aperture  function.  The  result  is  often  referred  to  as  the  presampling  NPS: 

Spre(u,v)  =  S(u,v)|absinc(au)sinc(bv)|  (17) 

where  S(u,v)  is  the  analog  NPS,  and  a  and  b  are  the  pixel  aperture  dimensions. 
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,  ’The  result.of  b)  is  that  the  presampling  NPS  is  convoluted  with  the  transform  of  the  sampling  comb  function: 

Sd(uIv)  =  Spre(uIv)  ®®  |||(u, v, (18) 

Ax  Ay 

where  Ax  and  Ay  are  the  sampling  intervals  in  the  x  and  y  dimensions,  and  ®  ®  denotes  two  dimensional  convolution. 

Note  that  this  comb  function  consists  of  an  infinite  array  of  delta  functions  separated  in  each  dimension  in  frequency 
space  by  two  times  the  respective  Nyquist  frequencies;  Un  =  ViAx,  vN  =  V6Ay.  S<jjg(u,v)  is  therefore  an  infinite  2-d  array  of 
replications  of  Spre(uPv),  with  the  origin  of  each  Spre  separated  by  2uN  in  the  u  direction  and  2vN  in  the  v  direction.  This 
replication  means  that  noise  beyond  the  Nyquist  frequency  overlaps  noise  below  the  Nyquist  frequency  of  adjacent 
replications.  This  overlap,  or  aliasing,  becomes  increasingly  severe  as  the  sampling  interval  is  increased,  that  is  as  Spre  is 
increasingly  undersampled.  Because  the  phases  of  the  Fourier  components  of  the  image  noise  are  random,  the 
presampling  NPS  cannot  be  directly  measured  using  fine  sampling  techniques  such  as  those  commonly  employed  to 
measure  the  presampling  MTF  [H.  Fujita,  D.  Tsai,  T.  Itoh,  K.  Doi,  J.  Morishita,  K.  Ueda,  and  A.  Ohtsuka,  “A  simple  method 
for  determining  the  modulation  transfer  function  in  digital  radiography”,  IEEE  Trans.  Med.  Imag.  11(1),  34-39  (1992)]. 

Using  digitized  screen-film  images,  Giger  et  al.  have  studied  the  effects  on  the  digital  NPS  of  the  sizes  of  the  integrating 
pixel  aperture  and  the  sampling  interval  [M.L.  Giger,  K.  Doi,  and  C.E.  Metz,  “Investigation  of  basic  imaging  properties  in 
digital  radiography”,  Med.  Phys.  1 1(6),  797-805, 1984.].  Their  results  show  that  aliasing  of  high  frequency  noise  (due,  for 
example,  to  the  high  frequency  components  of  x-ray  quantum  noise)  above  the  Nyquist  frequency  can  form  a  significant 
component  of  the  measured  digital  NPS.  Large  integrating  aperture  size  (i.e.  larger  than  the  sampling  interval)  tends  to 
reduce  the  effect  of  this  aliasing  by  smoothing  the  high  frequency  x-ray  quantum  noise,  at  the  expense  of  spatial 
resolution.  Aperture  sizes  less  than  the  sampling  distance  increase  the  overall  noise  level,  and  particularly  so  at  high 
frequencies.  Dobbins  has  recently  discussed  the  impact  of  undersampling  on  the  interpretation  of  the  MTF  of  digital 
systems  [J.T.  Dobbins,  “Effects  of  undersampling  on  the  proper  interpretation  of  modulation  transfer  function,  noise  power 
spectra,  and  noise  equivalent  quanta  of  digital  imaging  systems”,  Med.  Phys.  22(2),  171-181  (1995)]. 

The  result  of  c)  is  that  the  NPS  is  also  convoluted  with  the  square  of  the  transform  of  the  data  window,  w(x,y): 

Sd,w(u.v)  =  Sd(u,v)®®|W(4v)|2  (19) 

where  Sw(u,v)  is  the  measured  NPS,  and  the  frequency  window,  W(u,v)  is  the  Fourier  transform  of  the  data  window,  w(x,y) 
[N.C.  Geckinli  and  C.  Yavuz,  Discrete  Fourier  Transformation  and  Its  Application  to  Power  Spectra,  (Elsevier,  Amsterdam 
1983)  pg.  99]. 

If  the  data  series  is  simply  truncated  after  N  x  M  points,  corresponding  to  multiplication  of  the  infinite  data  sequence  by  a 
unity  height  box  of  dimensions  NAx  x  MAy,  then  the  NPS  is  convoluted  with  a  frequency  window  with  modulus: 

|W(u.«)|  =  NA»l^yl^uW.’,^a,->  (20) 

ir  NAxuMAyv 

A  variety  of  data  windows  have  been  studied  in  many  signal  processing  contexts,  with  the  optimum  choice  for  a  given 
situation  depending  on  the  characteristics  of  the  particular  spectrum  being  characterized.  Selection  of  a  data  window  is 
based  on  the  shape  and  width  of  the  resulting  frequency  window.  In  general,  there  is  a  compromise  between  the  width  of 
the  central  lobe  and  the  magnitude  and  extent  of  any  side  lobes.  A  broad  central  lobe  reduces  spectral  resolution,  and 
significant  side  lobes  produce  leakage,  wherein  noise  power  from  nearby  frequencies  are  mixed  in  the  convolution 
process.  In  the  case  of  power  spectra  exhibiting  strong  periodic  components,  as  is  the  case  in  many  digital  radiographic 
systems,  minimization  of  side  lobe  magnitude  is  often  necessary  for  distinguishing  small  amplitude  features  adjacent  to 
large  amplitude  peaks  [N.C.  Geckinli  and  C.  Yavuz,  Discrete  Fourier  Transformation  and  Its  Application  to  Power  Spectra, 
(Elsevier  Science  Publishing  Co.  Inc.,  Amsterdam,  1983)  Chap.  3.,  R.W.  Ramirez,  The  FFT:  Fundamentals  and  Concepts, 
(Tektronix,  Inc.,  Englewood  Cliffs,  NJ,  1985)  Chap.  7.]. 


D.  Experimental  Estimation  of  the  Noise  Power  Spectrum 

Because  noise  in  a  radiographic  image  arises  from  stochastic  (random)  processes,  any  realization  (image)  of  the  noise 
associated  with  those  processes  is  merely  a  sample  chosen  from  the  ensemble  of  all  possible  samples.  Thus  there  is  an 
inherent  statistical  uncertainty  in  the  representation  of  the  underlying  (true)  population  variance  by  that  sample.  In  this 
regard,  Fourier  decomposition  of  image  noise  (such  as  is  performed  in  determination  of  the  NPS)  is  different  from  the 
decomposition  of  signal  (such  as  is  performed  in  determination  of  the  modulation  transfer  function  (MTF)),  where 
amplitudes  and  phases  are  known  and  fixed.  The  variance  in  the  spectral  components  of  the  measured  NPS  is 
independent  of  the  number  of  data  points  (e.g.  the  number  of  pixels)  in  the  sample  used  to  calculate  it,  and  for  mean-zero, 
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.  ’  normally  distributed  noise,  results  in  a  coefficient  of  variation  (COV)  of  approximately  unity  for  each  spectral  estimate, 
S(u,v).  [R.F.  Wagner  and  J.M.  Sandrik,  “An  introduction  to  digital  noise  analysis”,  from  The  Physics  of  Medical  Imaging: 
Recording  System  Measurements  and  Techniques,  A.  Haus,  ed.  (AIP,  New  York,  1979),  Appendix  A].  The  variance 
(uncertainty)  in  the  estimate  can  be  reduced  by  sectioning  the  entire  NAx  x  MAy  noise  record  into  a  number  of  contiguous 
blocks,  and  averaging  the  spectra  obtained  from  each  (ensemble  averaging).  The  COV  of  the  final  estimate  decreases  as 
the  inverse  square  root  of  the  number  of  determinations.  However,  the  spectral  resolution  is  inversely  proportional  to  the 
number  of  data  points  used  to  generate  each  NPS.  Thus  for  a  fixed  amount  of  data,  there  is  a  tradeoff  between  frequency 
resolution  and  variance  in  the  final  NPS  estimate.  For  example,  if  N  data  points,  spaced  at  an  interval  Ax,  are  used  to 
calculate  a  1-d  NPS,  then  the  spectral  resolution  and  COV  are  respectively  [R.F.  Wagner  and  J.M.  Sandrik,  “An 
introduction  to  digital  noise  analysis",  from  The  Physics  of  Medical  Imaging:  Recording  System  Measurements  and 
Techniques,  A.  Haus,  ed.  (AIP,  New  York,  1979)]: 


COV  =  10 


If  these  N  points  are  divided  into  n  successive  sections,  then  the  frequency  resolution  and  COV  become: 


Af  = 


I  n 
f^Ax-NAx 


(21) 


(22) 


III.  METHODS 

A.  Data  Acquisition 

Both  mammographic  systems  utilize  molybdenum  target  x-ray  tubes,  with  molybdenum  or  rhodium  filtration.  Data  were 
acquired  using  Mo/Mo  spectra  and  a  variety  of  kVp  and  mAs  settings.  A  block  of  acrylic  3.8  cm  thick  and  sufficiently  large 
to  cover  the  entire  detector  surface  was  used  to  create  realistic  spectral  composition  and  scatter  conditions.  Detector 
entrance  exposure  was  measured  using  an  ion  chamber  placed  behind  the  acrylic  block,  with  corrections  applied  to 
account  for  slight  differences  between  locations  of  the  chamber  and  the  detector  entrance  surface,  relative  to  the  focal 
spot.  For  each  selected  kVp  and  mAs,  eight  images  were  obtained.  Results  given  here  are  for  data  taken  from  a  selected 
CCD  module  (i.e.  the  data  do  not  span  module  boundaries,  nor  are  data  from  different  modules  of  a  given  detector 
averaged).  The  data  blocks  were  divided  into  the  maximum  number  of  contiguous  sections  consistent  with  the  NPS 
estimation  technique  to  be  used  (e.g.  N  x  L  pixel  sections  for  a  scan  in  the  x-direction,  or  N  x  M  pixels  for  a  2-d  NPS 
calculation).  Spectra  calculated  from  each  section,  and  for  each  of  the  eight  images  were  then  averaged.  The  final 
ensemble  averaged  NPS  estimates  were  thus  averages  of  -  300  - 1200  spectra,  resulting  in  standard  errors  of  -3%-  6%. 


B.  Subtraction  of  Background  Trends 

Low  frequency  background  trends,  such  as  those  from  the  x-ray  source  heel  effect,  are  often  present  in  the  data.  Although 
such  trends  may  have  significant  Fourier  components  only  at  frequencies  below  the  lowest  measurable  frequency  in  the 
power  spectrum  (1/NAx),  leakage  into  low  frequency  spectral  components  can  artificially  inflate  the  low  frequency  NPS. 
For  this  reason,  noise  power  spectra  are  frequently  obtained  by  first  subtracting  one  uniform  exposure  image  from  another 
made  with  the  same  exposure,  and  dividing  the  NPS  resulting  from  the  subtracted  image  by  two.  The  factor  of  two 
accounts  for  the  fact  that  the  variance  in  the  subtracted  image  gets  approximately  equal  contributions  from  the 
uncorrelated  noise  in  the  two  original  images.  Alternately,  a  difference  image  may  be  obtained  by  subtracting  a  highly 
averaged  uniform  exposure  image  from  a  single  image.  The  random  noise  in  the  highly  averaged  image  is  negligible  if  a 
sufficient  number  of  images  are  used  to  obtain  it. 

The  image  subtraction  technique  has  the  drawback,  however,  of  also  removing  any  uncorrected  fixed  pattern  noise 
occurring  at  higher  frequencies.  Thus,  this  method  may  not  be  optimum  for  measuring  the  NPS  of  images  that  have 
already  been  uniformity  corrected  (flat  fielded),  because  noise  removed  in  the  subtraction  process  would  nevertheless  be 
present  in  clinical  images.  We  have  therefore  chosen  not  to  use  image  subtraction,  but  to  remove  low  frequency 
background  trends  in  other  ways.  Two  methods  have  been  evaluated:  subtraction  of  a  low  pass  filtered  version  of  the 
image,  and  subtraction  of  a  low  order  polynomial  fit  to  the  image.  When  using  the  first  method,  for  simplicity  a  rectangular 
(boxcar)  filter  was  used  to  filter  the  entire  2-d  data  set.  Subtraction  of  the  resulting  filtered  image  from  the  original  is 
equivalent  to  subtraction  of  a  running  local  mean  as  has  been  employed  for  one-dimensional  data  by  others  [M.L.  Giger, 
K.  Doi,  and  C.E.  Metz,  “Investigation  of  basic  imaging  properties  in  digital  radiography.  2.  Noise  Wiener  spectrum",  Med. 
Phys.  'l  1  (6),  797-805  (1984).].  Convolution  in  coordinate  space  with  a  rectangular  filter  is  equivalent  to  multiplication  in 
frequency  space  by  a  sine  function.  The  dimensions  of  the  square  filtration  kernel  were  chosen  sufficiently  large  that  the 
first  zero  crossing  of  the  resulting  sine  function  was  at  or  below  the  minimum  frequency  sampled  in  the  NPS.  In  the 
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.  ’  second  method  of  background  trend  subtraction,  a  2-dimensional  polynomial  was  first  fit  to  each  uniform  exposure  image, 
and  subtracted.  Both  first  and  second  degree  polynomials  were  evaluated. 

In  order  to  compare  the  above  two  methods  of  background  trend  subtraction,  a  synthetic  data  set  was  created.  The  data 
consist  of  a  DC  (average)  pixel  value,  a  normally  distributed  random  component,  a  periodic  component  (one  in  each 
dimension),  a  linear  (ramp)  component,  and  a  2™  order  component.  The  digital  value  of  pixel  l(i,j)  was  computed  as: 

l(i,j)  =  d1*[rand(i,j)]  +  d2(i,j)  (23) 

where  the  rand()  operator  generates  one  of  a  set  of  normally  distributed  random  numbers  with  average  zero  and  standard 
deviation  one,  dT  is  a  scaling  factor,  and 

d20.j)  =  aosin(iup+  8u)  +  b0sin(jvp+  6v)  +  a1i+b1j+a2i2  +b2j2  +C0  (24) 

where  up  and  vp  are  the  spatial  frequencies  of  the  periodic  components  in  the  x-  and  y-directions,  respectively,  and  8„  and 
5v  are  phase  factors.  The  constant  Co  sets  the  magnitude  of  the  DC  component,  and  the  constants  ao,  bo,  a1f  b1t  a2,  and  b2 
set  the  magnitudes  of  the  periodic,  linear,  and  quadratic  components  in  the  x-  and  y-directions,  respectively.  Figure  1  is  a 
grayscale  plot  of  a  portion  of  one  of  the  synthetic  data  sets. 

Figure  2  shows  a  semilog  plot  of  four  representative  1-d  power  spectra  generated  from  synthetic  data.  The  spectrum 
labeled  with  diamonds  is  that  of  a  data  set  with  no  background  trend  (a1  =  a2  =  bi  =  b2  =  0).  The  second  spectrum,  labeled 
with  squares,  is  from  the  same  data,  but  with  linear  and  quadratic  components  added.  The  third  and  fourth  spectra 
(triangles  and  crosses),  are  calculated  following  trend  removal  via  subtraction  of  a  least  squares  surface  fit  to  a  second 
order  polynomial  or  by  subtraction  of  a  boxcar  filtered  version  of  the  image,  respectively.  In  this  example,  the  boxcar  filter 
was  64  x  64  pixels,  the  same  size  as  the  sections  of  data  used  to  calculate  the  2-d  NPS.  The  spectrum  obtained  with  no 
trend  removal  demonstrates  substantial  leakage  of  low  frequency  noise  power  into  the  frequency  range  of  the  periodic 
component.  The  figure  shows  that  both  trend  removal  techniques  reduce  the  noise  power  sufficiently  at  low  frequency  to 
eliminate  the  effects  of  the  leakage.  The  quadratic  fit  technique  virtually  restores  the  original  NPS,  a  result  not  surprising 
given  the  quadratic  nature  of  the  trend. 

Figure  3  also  shows  the  effects  of  different  trend  removal  techniques,  this  time  using  flat  field  image  data  obtained  from 
one  of  the  mammography  detectors.  The  2-d  spectra  were  computed  from  the  average  of  648  128  x  128  pixel  sections. 
Thin  slices  were  extracted  along  the  v-axis.  Again,  slices  are  shown  from  2-d  spectra  following  no  trend  removal,  and 
following  subtraction  of  either  a  fitted  quadratic  surface  or  a  low-pass  filtered  data  set.  In  addition,  two  different  size  . 
convolution  kernels  for  low  pass  filtering  are  compared.  For  a  rectangular,  n  x  m  pixel  convolution  kernel  of  amplitude 
(nAxmAy)'1,  where  Ax  and  Ay  are  the  pixel  size  in  the  x-  and  y-dimensions,  the  trend-corrected  NPS  is  given  by: 

®d,w,b  =  Sd,w[l  -  sinc(nAxu)  sinc(mAyv)]  (25) 

In  order  that  only  noise  power  close  to  DC  is  removed,  one  would  like  the  term  sinc(nAxu)sinc(mAyv)  to  be  appreciable 
only  for  very  low  frequency.  Therefore,  the  quantities  nAx  and  mAy,  corresponding  to  the  reciprocal  of  the  first  zero 
crossings  of  the  sine  functions  in  the  two  dimensions,  must  be  sufficiently  large  compared  to  the  lowest  frequency  of 
interest  in  the  NPS.  In  Figure  3,  SdiWib  is  shown  for  n=m=100,  and  for  n=m=200.  Note  that  for  m=100,  the  first  zero 
crossing  of  sinc(mAyv)  is  at  a  slightly  higher  frequency  than  the  lowest  frequency  (v=1/(128Ay))  in  the  spectrum.  The 
effect  is  a  slight  reduction  in  the  noise  power  at  the  lowest  frequency,  relative  to  that  shown  for  m=200.  Because  of  this,  in 
general  we  have  found  the  best  results  by  making  the  size  of  the  smoothing  kernel  somewhat  larger  than  the  size  of  the 
blocks  used  to  calculate  the  individual  2-d  spectra. 

It  should  be  noted  that  2-d  quadratic  fitting  carries  a  significantly  higher  cost  in  processing  time  than  does  2-d  smoothing. 
This  is  true  even  if  the  fits  are  performed  on  the  relatively  smaller  blocks  individually,  rather  than  on  the  full  data  set.  For 
example,  on  a  Sun  Sparcstation  10  in  our  lab,  trend  removal  using  2-d  fits  to  8Q0  128  x  128  blocks  takes  approximately  4 
hours  of  processing  time,  whereas  the  same  amount  of  data  can  be  smoothed  in  several  minutes. 

C.  Units  and  Scaling 

The  noise  power  spectrum  is  more  precisely  termed  the  power  spectral  density,  that  is  the  noise  power  per  unit  frequency. 
[R.B.  Blackman  and  J.W.  Tukey,  The  Measurement  of  Power  Spectra,  (Dover  Publications,  Inc.,  New  York,  1958),  Chap. 
1].  The  quantity  S(u,v)  du  dv  is  the  contribution  to  the  variance  from  noise  components  with  spatial  frequencies  between  u 
and  u  +  du,  and  v  and  v  +  dv.  The  correct  scaling  of  the  NPS  can  be  verified  by  observing  that  the  volume  (area)  under 
the  2-d  (1-d)  NPS  equals  the  total  variance.  Thus  integration  over  frequency  of  spectra  estimated  using  mean-zero 
sections  from  uniform  illumination  images  should  numerically  equal  the  mean  square  deviations  from  their  average  of  the 
pixel  values  in  those  sections. 
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Noise  power  spectra  are  sometimes  presented  normalized  by  the  square  of  the  large  area  (average)  output  signal  [J.T. 
Dobbins  III,  D.L.  Ergun,  L.  Rutz,  D.A.  Hinshaw,  H.  Blume,  and  D.C.  Clark,  "DQE(f)  of  four  generations  of  computed 
radiography  acquisition  devices",  Med.  Phys.  22(10),  1581-1593  (1995),  R.J.  Jennings,  H.  Jafroude,  R.M.  Gagne,  T.R. 
Fewell,  P.W.  Quinn,  D.E.  Steller  Artz,  J.J.  Vucich,  M.T.  Freedman,  and  S.K.  Mun,  “Storage-phosphor-based  digital 
mammography  using  a  low-dose  x-ray  system  optimized  for  screen-film  mammography",  SPIE  2708,  220-232  (1996).]. 
This  normalization  has  the  effect  of  compensating  for  differences  in  gain  between  two  systems  being  compared,  and  of 
simplifying  the  expression  for  the  DQE.  With  this  normalization,  equation  (5)  becomes: 


DQE(u,v,<p) 


C(u,v)MTF2(u,v) 
§(u,v,<p)  SNR2 


(26) 


where  S(u,v,<p)  =  S(u,v,<p)/((pg)2  is  the  normalized  NPS. 


IV.  RESULTS 

A.  2-dimensional  NPS 

The  additive  noise  in  the  NPS  obtained  from  a  CCD-based  digital  mammography  detector  can  be  written  as: 

ai  =  a?+a2  +  Sp(u,v)  (27) 

where  ot  is  the  CCD  thermal  (Johnson)  noise,  <rr  is  the  read  noise,  and  Sp(u,v)  is  the  variance  due  to  frequency  dependent 
periodic  noise.  Periodic  noise  can  arise  from  capacitively  coupled  clock  pulses,  electromagnetic  interference  from  power 
supplies,  or  60  Hz  pickup  due  to  ground  loops.  Because  the  phase  (in  time)  of  such  periodic  interference  is  not  fixed  with 
respect  to  detector  readout,  its  phase  in  the  image  is  not  fixed  from  frame  to  frame,  and  it  is  therefore  not  removed  during 
the  uniformity  correction. 

2-d  spectra  obtained  from  both  developmental  systems  exhibit  strong  periodic  structures  over  all  exposure  ranges  tested. 
These  structures  take  the  form  of  sharp  peaks  in  the  spectra,  and  are  due  to  periodic  additive  noise.  Figure  4  shows 
grayscale  plots  of  the  2-d  spectra  of  one  digital  mammographic  system  at  four  exposure  levels.  Gridlines  are  spaced  at 
one  mm'1.  With  increasing  x-ray  exposure,  the  periodic  noise  peaks  are  decreasingly  apparent  as  the  grayscale  level  is 
changed  to  accommodate  the  higher  x-ray  quantum  noise  contribution  at  low  frequency. 


B.  Effect  of  Varying  Slit  Length  on  1-d  NPS 

We  have  studied  the  effect  of  varying  the  length  of  the  synthetic  scanning  slit  used  for  the  estimation  of  S(u,0)  or  S(0,v). 
Figure  5  shows  1-d  spectra  obtained  from  uniform  illumination  images  using  an  exposure  of  38  mR  at  the  detector 
surface.  Each  of  the  spectra  shown  in  the  plot  was  generated  by  stepping  a  synthetic  slit  with  x  and  y  dimensions  of  1  x  L 
pixels,  over  the  data  in  the  x-direction  in  steps  of  one  pixel.  The  longest  slit  length  shown  (L=512  pixels)  corresponds  to  a 
slit  approximately  20  mm  long.  At  each  step,  the  pixels  in  the  slit  were  averaged.  The  resulting  128  element  series  was 
then  Fourier  transformed  and  scaled  as  follows  to  give  the  measured  1-d  NPS.  Combining  equations  (14)  and  (15)  we 


have: 


S(kAu)=  Jim  NLAy^ 

NAx-*~  N 

NLAy — 


I  |'(nAx)e-2wi<kAu)(nAx) 
n=0 


|2 


(28) 


1  Nl-1 

where  l'(nAx)  =  —  £l(nAx,IAy)  is  the  average  of  the  pixel  values  in  the  synthetic  slit.  Note  that  from  equation  (16),  a 
Nl  i=o 

section  of  the  digital  2-d  NPS  along  the  u-axis,  calculated  from  a  N  x  M  section  of  the  image,  is: 

|2 


S(kAu,0)  =  lim  ^TTr 
NAx->~  N  M 
MAy-»~ 


N-1 

I  e 

n=0 


-27t{(kAu)(nAx)] |(nAx  mA  J 


(29) 


m=0 


For  the  special  case  when  the  slit  length,  NLAy,  is  equal  to  the  scan  length,  NAx,  equations  (28)  and  (29)  reduce  to  the 
same  expression: 


S(kAu,0)  = 


NAx->oo  N  N{_ 
N^Ay— 


N-1 
I  e 


-2jt((kAu)(nAx)]  r(nAx) 


n=0 


(30) 
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Thus,  for  example,  a  cut  along  the  u-axis  from  a  128  x  128  2-d  NPS  is  identical  to  the  1-d  NPS  obtained  by  scanning  a  1  x 
128  pixel  synthetic  slit  over  128  pixels  in  the  x-direction.  In  all  cases,  we  verified  that  the  1-d  NPS  obtained  from  the 
scanned  slit  method  with  the  slit  length  equal  to  the  size  of  the  data  block  used  to  generate  the  2-d  NPS,  was  identically 
equal  to  the  section  from  the  2-d  spectrum  along  the  u-  or  v-axis. 

As  shown  in  Figure  5,  the  very  low  frequency  values  of  the  1-d  NPS  did  not  converge  to  a  plateau  value  for  slit  lengths  of 
up  to  512  pixels,  which  is  twice  the  scanned  distance.  This  behavior  is  similar  to  that  observed  by  Sandrik  and  Wagner  for 
the  Kodak  Hi-Plus/XRP  screen-film  system.  [J.M.  Sandrik  and  R.F.  Wagner,  “Radiographic  screen-film  noise  power 
spectrum:  variation  with  microdensitometer  slit  length”,  Appl.  Opt.  20(16),  2795-2798  (1981).].  In  that  study,  measured 
NPS  values  at  0.39  cycles/mm  continued  to  increase  with  increasing  slit  length  up  to  L  =  5.9  mm  for  scan  a  distance  of  2.5 
mm,  while  convergence  to  a  plateau  value  was  observed  for  other  screen-film  systems  at  similar  frequencies. 

C.  Slices  from  the  2-d  NPS 

Notably  in  the  context  of  the  evaluation  of  storage  phosphor  systems,  others  have  obtained  1-d  spectral  estimates  by 
extracting  thick  cuts  from  the  2-d  NPS.  The  cuts  are  made  parallel  to  the  primary  axis  of  interest,  but  do  not  include  the 
axis,  so  as  to  avoid  low  frequency  trending  effects.  We  have  evaluated  this  approach  for  analysis  of  the  spectra  of  the 
digital  mammography  systems.  Figure  6  is  a  plot  comparing  a  thin  cut  (a  single  frequency  bin  wide)  along  the  primary  (v) 
axis  of  the  2-d  NPS  with  a  thick  cut  generated  by  summing  the  two  four-bin-wide  sections  immediately  on  either  side  of, 
and  adjacent  to,  the  v-axis.  The  thick  slices  adequately  represent  the  on-axis  slice  except  at  low  frequency,  where  the 
NPS  decreases  rapidly  with  increasing  u.  Of  greater  significance  is  the  omission  of  the  periodic  noise  structures  near 
spatial  frequencies  of  3  and  6  cycles/mm.  These  represent  noise  power  from  fluctuations  in  the  y-direction  (left  to  right 
direction  of  the  breast  support)  that  are  constant  in  the  anterior-chest  direction.  Such  noise  may  appear  as  horizontal 
stripes  in  the  mammographic  image. 

As  would  be  expected,  thicker  slices  are  poorer  approximations  to  the  on-axis  cuts,  as  the  falloff  in  the  2-d  NPS  in  the 
dimension  orthogonal  to  the  slices  causes  the  average  spectral  values  to  be  reduced.  This  produces  an  effect  similar  to 
that  resulting  from  utilization  of  a  scanning  slit  of  insufficient  length. 

D.  NPS  of  the  Developmental  Mammographic  Systems 

The  NPS  of  uniform  illumination  images  from  both  digital  mammographic  systems  was  measured  over  a  range  of  detector 
entrance  surface  exposures  from  9  mR  to  78  mR.  For  each  exposure,  2-d  noise  power  spectra  estimates  were  made  as 
described  in  Section  III. A.  Figure  7a)  and  7b)  are  semilog  plots  of  cuts  along  the  u-  and  v-axes  of  both  systems  at  four 
exposure  levels.  The  spectra  have  been  normalized  by  the  respective  large  area  output  signal  squared.  Immediately 
apparent  are  the  very  different  curvatures  of  the  spectra  obtained  from  the  two  systems,  with  the  spectra  of  System  1 
leveling  off  at  high  frequency,  while  those  of  System  2  continue  decreasing  out  to  the  Nyquist  frequency.  At  all  exposures 
tested,  the  normalized  noise  power  of  System  2  is  less  than  Iff6  at  Nyquist.  Also  apparent  from  Figure  7  is  the  fact  that 
the  noise  behavior  of  System  1  is  different  in  the  x-  and  y-dimensions,  with  the  y-dimension  exhibiting  significantly  greater 
noise  at  high  frequency.  In  both  systems,  the  y  direction  corresponds  to  the  parallel  transfer  direction  of  the  CCDs. 


E.  Comparison  with  Pixel  Variance 

Each  of  the  2-d  spectra  presented  in  this  study  was  integrated  over  spatial  frequency  in  order  to  compare  the  resulting 
value  with  the  variance  calculated  from  the  pixel-to-pixel  fluctuations.  The  integrals  were  calculated  over  the  range  from  - 
fN  to  +fN,  where  fN  is  the  Nyquist  frequency.  In  all  cases,  differences  between  the  integrals  and  the  directly  calculated  total 
variance  were  less  than  0.5%. 

Figure  8  shows  the  total  noise  of  the  two  digital  mammography  detectors  over  the  range  of  exposures  from  9  mR  to  78 
mR.  The  noise  was  determined  from  the  square  roots  of  the  integrated  volumes  under  the  2-d  spectra  whose  axial  slices 
are  shown  in  Figure  7.  The  figure  plots  the  log  of  the  noise  versus  the  log  of  the  exposure.  For  both  systems,  straight  lines 
have  been  fit  to  the  experimentally  determined  data  points,  and  the  calculated  slopes  of  the  lines  are  shown  in  the  figure. 
In  a  system  for  which  x-ray  quantum  noise  is  the  dominant  noise  source,  the  slope  should  be  approximately  0.5. 


V.  DISCUSSION 

Similarly  to  other  types  of  radiographic  systems,  digital  detectors  for  mammography  can  usefully  be  characterized  through 
Fourier  analysis  of  their  noise  characteristics.  However,  investigators  preparing  to  estimate  noise  power  spectra  of  such 
systems  are  confronted  with  a  bewildering  array  of  procedural  considerations,  ignorance  of  which  can  lead  to  misleading 
results.  Because  significant  digital  processing  power  is  now  available  to  the  majority  of  medical  physicists,  limitations 
formerly  imposed  by  the  tradeoff  between  spectral  resolution  and  uncertainty  in  the  spectral  estimates  are  much  less 
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,  stringent  tfian  even  a  decade  ago.  Thus  ensemble  averaging  over  a  large  number  of  uniform  exposure  images  is  feasible, 
thereby  permitting  fairly  subtle  features  of  the  NPS  to  be  resolved.  The  necessary  resolution  is  determined  by  the 
proximity  in  Fourier  space  of  relevant  features  and  should  be  determined  on  a  system-by-system  basis.  For  the  digital 
systems  characterized  in  this  paper,  128  x  128  2-d  spectra,  yielding  64  frequency  bins  in  the  single-sided  1-d  spectra, 
were  sufficient  to  resolve  the  features  of  interest.  Occasionally,  we  found  it  useful  to  recalculate  particular  spectra  at 
higher  resolutions.  This  may  be  useful,  for  example,  for  assurance  that  low  frequency  structures  are  in  fact  independent  of 
background  trend  noise  power  residing  in  the  lowest  frequency  bin  (the  range  between  zero  and  1/NAx,  where  NAx  is  the 
length  of  the  data  series). 

The  noise  power  in  images  from  both  the  digital  mammographic  systems  has  markedly  non-isoplanatic  properties, 
particularly  at  low  exposure  where  frequency-dependent  system  noise  is  most  predominant.  Because  the  inherent  nature 
of  the  CCD  readout  process  defines  orthogonal  serial  and  parallel  directions,  image  noise  generated  during  the  readout 
process  is  unlikely  to  fall  along  the  lul  =  Ivl  lines  of  the  2-d  NPS.  Such  noise  can  be  due  to  digital  clock  pulses  that  are 
radiatively  or  capacitively  coupled  to  the  analog  data  signals.  Furthermore,  the  majority  of  this  noise  power  resides  away 
from  the  principle  spatial  frequency  axes,  and  thus  does  not  appear  on  the  1-d  NPS,  whether  it  is  obtained  from  an 
extracted  slice  of  the  2-d  NPS,  or  from  a  scanned  slit.  It  is  therefore  important  that  the  full  2-d  NPS  be  estimated,  even  if 
1-d  spectra  are  subsequently  used  for  SNR  calculations. 

The  scanned  slit  method  for  obtaining  1-d  spectra  yielded  results  identical  to  those  obtained  from  thin  cuts  through  the  2-d 
spectra,  except  for  the  susceptibility  of  the  former  approach  to  slit-length-dependent  variability  in  the  low  frequency 
components.  The  absence  of  convergence  of  the  low  frequency  values  of  S(u,0)  reflects  the  peaked  nature  in  the  v- 
direction  of  the  2-d  NPS,  S(u,v),  near  v=0.  Since  the  value  of  S^u)  obtained  from  scanning  is  scaled  by  L  to  get  S(u,0) 
(c.f.  equation  (14)),  only  for  values  of  u  where  the  addition  to  the  integral  of  equation  (13)  is  proportional  to  L  will  the  result 
be  insensitive  to  changes  in  L.  Thus,  approximating  the  function  sinc?(Lv)  as  a  triangle  of  base  width  2/L  and  unity  height 
centered  at  v=0,  S(u,0)  obtained  using  a  slit  of  length  Li  is  approximately  equal  to  that  obtained  with  slit  length  L2  only  if 
S(u,v)  is  approximately  constant  in  the  v-direction  between  v=2/L1  and  v=2/L2.  Note  that,  except  for  u  =  0,  where  S(0,v)  is 
symmetric  about  the  u-axis,  this  requirement  must  also  be  met  between  v=-2/Li  and  v=-2/L2.  Dainty  and  Shaw  have 
recommended  that  the  slit  length  be  at  least  equal  to  the  inverse  of  the  frequency  resolution  (that  is,  at  least  equal  to  the 
scan  distance).  Our  results  indicate  that  for  digital  systems,  a  length  significantly  greater  than  the  scan  distance  may  be 
required  for  asymptotic  convergence  of  the  low  frequency  NPS  estimates.  In  practice,  the  most  direct  method  of 
determining  sufficiently  large  L  is  probably  inspection  of  cross-sectional  cuts  through  the  2-d  NPS,  obtained  at  as  high  a 
frequency  resolution  as  is  practicable,  bearing  in  mind  the  guidelines  of  the  preceding  sentences. 

For  the  noise  data  studied  here,  subtraction  of  a  duplicate  data  set  that  has  been  convoluted  with  a  two-dimensional  rect 
function  kernel  is  an  effective  and  efficient  method  for  the  removal  of  low  frequency  trending  effects.  The  effectiveness  of 
the  technique  is  indicated  by  the  similarity  between  thick  slices  extracted  from  the  2-d  NPS  near  the  central  axes,  and  on- 
axis  slices  (c.f.  Figure  6).  Because  the  2-d  NPS  can  contain  on-axis  features  not  present  in  adjacent  regions  of  frequency 
space,  1-d  spectra  obtained  from  thick  slices  parallel  to  the  primary  axis  of  interest  may  not  be  representative  of  the  on- 
axis  NPS.  For  this  reason,  thin  on-axis  cuts  are  recommended. 


VI.  CONCLUSIONS 

Fourier  analysis  has  been  used  to  characterize  the  noise  properties  of  two  developmental  detectors  for  digital 
mammography.  The  noise  power  spectra  of  images  from  the  two  systems  exhibit  several  common  features.  They  contain 
features  that  are  strongly  dimension-specific,  and  are  often  localized  to  relatively  small  portions  of  frequency  space.  They 
also  contain  predominant  peaks  indicative  of  periodic  additive  noise,  both  on  and  off  the  primary  frequency  axes.  On  the 
other  hand,  the  noise  properties  of  the  two  mammographic  detectors  differ  significantly  in  some  aspects.  Over  a  range  of 
input  exposures,  the  spectra  of  System  1  are  characterized  by  nearly  constant  noise  power  for  frequencies  greater  than 
about  8  cycles/mm,  while  those  of  System  2  exhibit  increasingly  sharp  fall  off  with  increasing  frequency.  Cuts  along  the 
two  primary  axes  of  the  2-d  NPS  of  System  1  exhibit  very  different  shapes,  while  those  of  System  2  are  nearly  identical 
along  the  two  axes.  At  all  frequencies,  the  relative  noise  of  System  2  decreases  monotonically  with  increasing  x-ray 
exposure,  while  that  of  System  1  reaches  a  minimum  then  begins  to  increase  again.  The  overall  noise  of  System  one 
increases  approximately  as  the  square  root  of  the  number  of  incident  x-ray  quanta,  while  that  of  System  2  increases  more 
slowly. 

Because  of  the  complex  nature  of  the  frequency  distribution  of  the  noise  power  in  images  from  digital  mammographic 
systems,  analysis  of  the  complete  2-d  spectra,  rather  than  only  the  1-d  spectra,  is  important.  Straightforward  removal  of 
background  trending  artifacts  is  possible  via  a  simple  filtering  process.  Cuts  obtained  along  the  primary  axes  of  the  2-d 
spectra  can  provide  valid  estimates  of  the  1-d  NPS  without  the  use  of  synthetic  scanned  slit  techniques. 
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Figure  1 

Grayscale  image  of  a  synthetic  data  set  used  to  evaluate  methods  of 
background  trend  removal.  For  ease  of  visualization,  the  data  set  shown 
contains  only  sinusoidal,  random,  and  DC  components  (no  linear  or  quadratic 
component). 
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Background  Trend  Removal 


Figure  2 

Comparison  of  background  trend  removal  techniques,  using  synthesized  data  set.  The  figure 
shows  the  NPS  calculated  from  a)  data  with  no  background  trend  (diamonds),  b)  the  same  data 
as  a),  but  with  linear  and  quadratic  trends  added  (squares),  c)  the  same  data  as  b),  but  after 
subtraction  of  a  2-dimensional  second  order  polynomial  fit  to  the  data  (triangles),  and  d)  the 
same  as  b),  but  after  subtraction  of  a  low-pass  filtered  version  of  the  data  (crosses). 
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Trend  Subtraction 


Spatial  Frequency  (mm*1) 


Figure  3 

Comparison  of  background  trend  removal  techniques,  applied  to  images  obtained  with  one  of  the 
digital  mammography  systems,  obtained  using  nominally  uniform  illumination.  The  graph  shows 
four  1-d  power  spectra,  S(0,v),  obtained  from  thin  cuts  from  four  different  2-d  spectra.  The  2-d 
spectra  were  calculated  following  a)  no  trend  removal  (none),  b)  subtraction  of  a  fitted  quadratic 
surface  (fit),  c)  subtraction  of  a  low  pass  filtered  version  of  the  data  using  a  100  x  100  pixel 
convolution  filter,  and  d)  same  as  c),  but  using  a  200  x  200  pixel  filter. 
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Figure  4a 

2-d  NPS  at  9  mR  exposure  to  the  detector 
surface.  Gridlines  are  spaced  by  1  mm'1. 


Figure  4c 

2-d  NPS  at  38  mR  exposure  to  the  detector 
surface.  Gridlines  are  spaced  by  1  mm'1. 


Figure  4b 

2-d  NPS  at  18  mR  exposure  to  the  detector 
surface.  Gridlines  are  spaced  by  1  mm'1. 
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Figure  4d 

2-d  NPS  at  76  mR  exposure  to  the  detector 
surface.  Gridlines  are  spaced  by  1  mm'1. 
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Scanned  Slit,  38  mR 


Figure  5 

1-d  noise  power  spectra  obtained  from  the  same  data  set,  but  with  various  scanning  slit  lengths. 
Spectral  components  above  approximately  0.6  cycles  per  millimeter  are  insensitive  to  changes  in 
slit  length  for  slit  lengths  of  64  pixels  or  greater.  However, 
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Figure  6 

Slices  from  2-d  NPS  of  digital  mammographic  system.  The  solid  line 
is  a  slice  one  bin  wide,  along  the  v-axis.  The  dotted  line  is  the  sum 
of  two  4-bin  sections  taken  from  either  side  of  the  v-axis. 

Differences  between  the  two  occur  primarily  at  low  frequency  and 
because  of  on-axis  periodic  noise  components. 


MEDICAL  PHYSICS 


Mark  B.  Williams,  Ph.D. 


v> 

(L 

Z 

1 


System  2 


Figure  7 

Noise  power  spectra  of  the  two  developmental  digital  mammography 
systems  studied.  Spectra  are  shown  at  four  exposures  for  each  system. 
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